Vapor— liquid coexistence of the Stockmayer fluid in nonuniform external fields 
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We investigate the structure and phase behavior of the Stockmayer fluid in the presence of nonuni- 
form electric fields using molecular simulation. We find that an initially homogeneous vapor phase 
undergoes a local phase separation in a nonuniform field due to the combined effect of the field 
gradient and the fluid vapor-liquid equilibrium. This results in a high density fluid condensing in 
the strong field region. The system polarization exhibits a strong field dependence due to the fluid 
condensation. 



I. INTRODUCTION 

The application of external fields for manipulation of 
the structure and phase behavior of dipolar fluids has at- 
tracted growing interest in recent years jU [5] . Typically, 
both theoretical and experimental studies consider uni- 
form applied fields. However, in complex systems such as 
microfluidic devices field gradients occur naturally. Mo- 
tivated by experimental work on the demixing of binary 
mixtures in field gradients [3J, we have previously stud- 
ied theoretically the application of nonuniform electric 
fields to pure fluids [4] and simple mixtures [5j [6] . An- 
other experimental realization of the ability of field gra- 
dients to promote phase transitions is the field induced 
crystallization of colloidal suspensions promoted by di- 
electrophoretic forces [THS]. 

In particular, we examined the effect of nonuniform 
fields on the vapor-liquid coexistence [4] by combining 
the simple van der Waals mean field theory with On- 
sager's theory of dielectrics to investigate polar and non- 
polar fluids. Our main finding was that above a critical 
field, in situations where the fluid is unperturbed by a 
uniform field, a nucleation of a gas bubble from the liquid 
phase or a liquid droplet from the vapor phase is induced 
by a nonuniform field. This phase separation transition 
is promoted by the dielectrophoretic force which favors a 
higher permittivity (density) fluid in the region of strong 
field. The resulting modification in the fluid phase dia- 
gram is considerably larger compared to uniform fields. 

A system which can be considered as an idealized man- 
ifestation of nonuniform fields is a grand canonical en- 
semble where a uniform field, E = const., is applied 
within the system volume but not in the material reser- 
voir, where E = 0. An example is a slit pore in equi- 
librium with a bulk fluid. Clearly, in real systems there 
exists an interfacial region at the pore edges where the 
field is nonuniform. The field effect in this type of system 
has been studied for one component fluids [TUl E] where 
it was found to gradually increase the fluid density within 
the pore moderately. A stronger field effect was found by 
Brunet et al. P21 H3] who studied mixtures. They ob- 
served that if the mixture has a demixing instability and 



one of the components is dipolar, the coupling to the ex- 
ternal field leads to a pore filling transition which allows 
a sensitive control of the pore composition. 

The goal of this paper is to study the structure of 
a dipolar vapor confined in a slit pore and exposed to 
nonuniform electric fields in the canonical ensemble via 
molecular simulation. We compare our results with mean 
field theory and discuss the consequences of the full de- 
scription of the vapor-liquid interface in finite systems. 
The paper is organized as follows: Section [H] discusses 



the simulation methods and model system. Section HI A 
compares results for the fluid in a uniform field with pre- 
vious studies. The results for a nonuniform field are com- 



piled in Section MB Conclusions are given in Section 

El 



II. SIMULATION METHODS 

We consider N spherical dipolar particles with a diam- 
eter a and a permanent dipole moment \x. The particles 
interact via the Stockmayer pair potential: 
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where r%j = J", — rj stands for the displacement vector 
of particles i and j, e denotes the Lennard- Jones (LJ) 
interaction parameter and £o is the vacuum permittivity. 
In the following we use reduced units; length: r* = r/a, 
temperature: T* — ksT/e, density: p* = pa 3 , dipole 
moment: p* = pj 'y 1 AixEQea 3 and external field: E* — 
E\J 'Aire oea 3 ■ Here, ks is the Boltzmann constant and we 
take er = e = 47T£o = 1. For brevity we omit the asterisk 
superscript henceforth. 

The phase diagram of bulk and confined Stockmayer 
fluids in a uniform field was determined using a Gibbs En- 
semble - Hybrid Monte Carlo scheme (GE-HMC). The 
classical Gibbs Ensemble Monte Carlo (GEMC) simula- 
tion [14] offers a simple method to determine the vapor- 
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liquid equilibrium densities and pressures with a single 
simulation run. This is achieved by equating the chemi- 
cal potentials and pressures of two simulation boxes using 
appropriate Monte Carlo particle and volume exchange 
moves, respectively [15] ■ The third type of moves per- 
formed are particle translations, rotations and other con- 
formational changes of single particles within the two 
boxes. In the GE-HMC variation, single particle moves 
are replaced by a collective Hybrid Monte Carlo (HMC) 
move [TBI HZ] • A single HMC cycle consists of three steps: 
first, particles of the current configuration, o, are assigned 
new momenta and angular velocities by sampling a Gaus- 
sian distribution corresponding to the desired temper- 
ature. Second, the new configuration, n, is generated 
from a short MD trajectory in the microcanonial ensem- 
ble. Lastly, the new configuration is accepted/rejected 
according to the Metropolis criterion: 



min(l,exp(-A J ff/T)) 



(2) 



where AH = H(n) — H(o) is the resulting change in the 
system Hamiltonian. Detailed balance is satisfied if the 
integration algorithm used for the MD trajectory is time- 
reversible and area-preserving |17j . which is fulfilled by 
a simple velocity- Verlet integrator. All MD trajectories 
were produced with the ESPResSo package [18] . The col- 
lective HMC moves allow to efficiently sample the high 
density liquid phase and complex molecular configura- 
tions Q2]. 

GE-HMC simulation cycles were conducted with 512 
Stockmayer particles. The total simulation consisted of 
10 4 cycles and observable sampling was done after 2500 
equilibration cycles. A single cycle was composed of 100 
MC moves where the probability of the move type was as 
follows: 0.8 for a particle exchange move, 0.15 for a HMC 
move and 0.05 for a volume exchange move. The number 
of time steps in a HMC move was 10. Both the time step 
and the attempted volume change were adjusted during 
equilibration such that approximately 50% of the moves 
were accepted. 

Simulations in a nonuniform field where performed in 
a simplified model system, see Fig. [T] Consider the fluid 
confined in a wedge condenser made up from two flat 
electrodes with a potential difference V across them and 
an angle j3 between them. In this geometry the field in 
the azimuthal direction is perpendicular to the field gra- 
dient in the radial direction. This simplifies the solution 
of Gauss' law since it implies Ve-E = 00]. The resulting 
electric field is: 



E(r) 



(3) 



where R\ is the radius of the inner condenser wall, Ri+r 
the radial distance from the imaginary meeting point of 
the electrodes and 9 is the azimuthal angle. We focus 
on a small angular section 69 far from the electrodes of 
the capacitor and therefore rewrite Eq. (|3| in terms of 



Cartesian coordinates: 
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where E = V/(/3Ri) is the maximal field at z = and 
Aq = 1/R\ is a constant characterizing the length scale 
of the field gradient; for Aq = the field is uniform. 

In this paper we study the Stockmayer fluid under the 
external field given in Eq. Q. The contribution of this 
field to the potential energy of a single particle is given 
by 



ur l = -Hi-Ei = - 



1 + A z t 



(5) 



where Ei is the field at the site of particle i which is 
a function of the particle coordinate Zi. Via the poten- 
tial energy we derive the additional force and torque on 
particle i due to the field, i.e, 
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Note that the force contribution, Eq. ([6]), vanishes in 
the case of a uniform field. Moreover, Eq. (JgJ) indicates 
that particles are drawn to the strong field region and 
feel a stronger force when aligned with the field. This 
may be thought of as the microscopic origin of the di- 
electrophoretic force. 

MD simulations of the Stockmayer fluid in a nonuni- 
form field were performed using the suitably modified 
ESPREesSo package [TS]. During the simulation the LJ 
potential was cut off at r c — 3 and the long range dipolar 
potential was evaluated using the dipolar P 3 M algorithm 
[20) with metallic boundary conditions. 




FIG. 1. The model condenser: a wedge made of two flat elec- 
trodes with a potential difference V across them and an angle 
P between them. Ri (R2) is the distance of the inner (outer) 
insulating wall from the imaginary meeting point of the elec- 
trodes, r is the radial distance from the inner wall (thick 
line). The resulting electric field E is along the 6 direction. 
In the simulations we approximate a small angular segment 
58 as a slit. 
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When a nonuniform field is applied to the fluid its 
translational invariance in the direction of the field gradi- 
ent is broken and therefore periodic boundary conditions 
(PBC) can not be used in this direction. However, the 
implementation of the P 3 M method requires that we em- 
ploy PBC in all directions. Therefore, in the simulations 
we model the condenser as an infinite slab with the two 
confining walls placed at z = and z = D < L, where L 
is the cubic simulation box length. The unwanted dipo- 
lar interactions between slabs replicated along the z di- 
rection are corrected using the dipolar layer correction 
method of Ref. [21 . This allows us to use a small gap of 
empty space in the simulation box. In order to isolate the 
field effect, we use for the fluid- wall interaction a purely 
repulsive LJ potential shifted and cut off at r c — 2 1 / 6 
(WCA potential). 

Simulations of the dipolar fluid in the slab where ini- 
tialized from random particle configurations. In the sim- 
ulations we employ a Langevin thermostat. A time step 
of At = 0.004 was used in all simulations. The simula- 
tions equilibration period varied from 1 x 10 5 — 6 x 10 s 
time steps, depending on the field strength. The struc- 
tural and dielectric properties of the system were then 
sampled every 200 time steps for at least 4 x 10 5 time 
steps. Time averaged quantities sampled are denoted by 
(...). 



III. RESULTS AND DISCUSSION 

A. Phase behavior of the Stockmayer fluid in a 
uniform field 

We first tested the applicability of our GE-HMC sim- 
ulation by comparing our results to available data on the 
Stockmayer fluid with p — 2. Here, the standard long 
range correction is applied to the L J interaction [15] . The 
resulting coexistence curve is shown in the inset of Fig. [2] 
We obtain a critical temperature T c — 2.05 and density 
p c = 0.301, in good agreement with the results of Van 
Leeuwen el al. [22] (T c = 2.06, p c = 0.289) and Kiyohara 
et al. [H] (T c = 2.05, p c = 0.306). The small differences 
in critical parameters are probably due to fact that un- 
like the works above, the LJ potential in this study is cut 
off at a fixed radius. 

Henceforth, we will focus in this work on the Stock- 
mayer fluid with p = 1.5. Vapor-liquid coexistence 
curves for such a bulk Stockmayer fluid with and with- 
out an external field are shown in Fig. [2] Here, no long 
range correction is applied to the LJ interaction since we 
intend to compare these results with those obtained in 
the slab geometry. In the absence of an external field we 
find T c — 1.59 and p c = 0.305. In accord with previous 
studies, we find when a uniform field is applied the unsta- 
ble region in the phase plane is increased [24H28] . This is 
due to the increased dipole-dipole interaction and corre- 
lation as the dipoles get aligned by the field [Ml HE] ■ In 
particular, simulations of the Stockmayer fluid in an ex- 
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FIG. 2. Vapor-liquid coexistence of the bulk Stockmayer fluid 
with p = 1.5. Squares: in the absence of an external field. 
Circles: with a uniform field Eo = 4. Inset: in the absence of 
an external field for p = 2. Dashed curves are fits to the law 
of rectilinear diameters using the Ising exponent f3 — 0.326. 
Critical points are marked by hollow symbols. 



ternal field found reasonable agreement with the Landau 
mean field theory [5S] for the field effect on the critical 
temperature [35]. For Eq = 4 we find that the critical 
temperature increases to T c = 1.83 and p c — 0.301. 

The value of p = 1.5 for the dipole moment was chosen 
since it is suited for description of both molecular fluids 
[5U] and dipolar colloidal suspension alike [2 . Further- 
more, we will from now on set T = 1.6 < T c . Hence, the 
value of dipolar coupling constant, A = p 2 /T, is A = 1.41. 
This means that the dipolar interaction at a distance a 
is comparable to both the thermal and LJ interactions. 
Note that a A value close to unity below the critical tem- 
perature can only be realized for relatively small values 
of /z [S]. 

The effect of confinement on the coexistence curve of 
the Stockmayer fluid is that of suppression of the unsta- 
ble region of the phase diagram. Studies conducted so far 
focussed on narrow slabs where this effect is large |32j . 
The finite-size effects for wider slabs where not studied 
since they are quite small and also computationally pro- 
hibitively expensive. 

An estimate of finite-size effects is provided in the inset 
of Fig. [3j which shows the density profiles, (p(z)), for 
N = 8192 confined particles at T = 1.6 and under an 
uniform field E — 4x, parallel to the slab walls. For an 
average fluid density p = 0.3 (slab width of D w 25), close 
to the critical density, the density profile exhibits vapor- 
liquid coexistence. The dashed horizontal line in the inset 
corresponds to the bulk liquid density of pi — 0.629 which 
is only slightly higher than the average density of « 0.59 
for the liquid in the slab. In contrast, slowly expanding 
the slab such that a fluid average density of p = 0.05 
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FIG. 3. Density profiles (p(z)) in a nonuniform field for a 
Stockmayer fluid with an average density p = 0.05 and a 
temperature T = 1.6. For the field we used a magnitude 
of E = 4 in Eq. Q. The solid curve shows results for 
N = 8192 particles using Ao = 0.1 in Eq. Q. For the dash- 
dot curve N = 16384 and A = 0.0794; for the dashed curve 
N = 48000 and Aq = 0.0555. Dotted curve: same as the 
dashed curve but replacing the LJ part of the potential by 
a WCA potential. Inset: density profiles in a uniform field 
E — Ax at T = 1.6. Two average densities, p = 0.3 and 0.05 
are presented corresponding to slab widths of D « 25 and 
D » 46, respectively. 



is finally obtained (slab width of L m 46) results in a 
homogeneous vapor phase, see the dash-dot curve in the 
inset of Fig. [3] This is expected since p = 0.05 is smaller 
than the vapor phase density of the bulk fluid p v — 0.053. 



B. The Stockmayer fluid in a nonuniform field 

The situation is markedly different when a nonuniform 
field is applied. The dashed curve in Fig. [3] gives the 
density profile for N = 8192 particles with an average 
density p = 0.05 under the nonuniform field given by Eq. 
Q with E = 4 and A = 0.1. This profile shows the 
condensation of a liquid-like layer from the homogeneous 
vapor phase in the strong field region. The high density 
layer of width « 2— 3cr is followed by a sharp interface and 
then a distinct vapor phase. The width of the liquid-like 
layer grows as the fraction of energetically costly inter- 
face molecules is reduced in larger systems. This effect is 
shown in the dash-dot and solid curves in Fig. [3]where we 
increase the number of particles while keeping the aver- 
age density the same. Here, since we scale the simulation 
box to keep the average density constant we also adjust 
the field in Eq. Q through A such that E(z = D) is 
kept constant. In addition, the liquid-like layer density 
also increases due to the decreased energetic cost of the 
interface. 
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FIG. 4. Density profiles (p(z)) in a nonuniform field of differ- 
ent magnitudes. The fluid maximal density is in the strong 
field region close to the wall (circle markers) and it increases 
with the field magnitude. Inset: density profiles for a nonuni- 
form field perpendicular to the the slab walls (see text). For 
all curves Ao = 0.1, the fluid average density is p — 0.05 and 
the temperature is T — 1.6. 



We explain this condensation by the fact that the 
nonuniform field is large only in the vicinity of z = 0. 
Hence, particles are drawn to this region where they gain 
energetically both by aligning in the stronger field and 
also from the LJ interaction which compensates for the 
loss in entropy. The attractive short range part in the 
interaction is important for the formation of a dense liq- 
uid layer. To show this we also performed a simulation 
where we replace the LJ part of the interaction by the 
purely repulsive WCA potential. The result for the dipo- 
lar WCA fluid is shown in the dotted curve of [3| clearly, 
only a moderate increase in the fluid density occurs and 
it follows the gradual decay of the field. 

We use the Stockmayer potential parameters for water 
|30j . which gives p — 1.56 for the fluid and the field mag- 
nitude of Eq = 4 corresponds to a maximal local field of 
4.1 V/nm. At least for water confined at the molecular 
scale such a field is not unusual [10 . Nonetheless, this 
field magnitude is still 5-10 times larger than the fields 
required to induce condensation in the mean field treat- 
ment [3]. The high field is a consequence of the costly 
interfacial region in the small system we simulate and 
is expected to be reduced in the thermodynamic limit 
N,D -> oo. 

Fig. [IJshows the density profiles obtained as a function 
of the nonuniform field magnitude, E Q . As Eq increases 
the condensation occurs rapidly starting at E > 2.5 al- 
beit gradually as as our system is finite. The inset of Fig. 
[4] gives a comparison between the field of Eq. Q and a 
nonuniform field of the same functional form and mag- 
nitude but perpendicular to the slab walls. This type of 



5 




(a) 30 



20 



10 



l(b) 



10 15 20 

z 



FIG. 5. Snapshots of a thin cross section through the y axis, 
focusing on the strong field region. Arrows indicate the dipole 
moment's x and z components scaled by a factor of 2. In (a) 
the electric field is parallel to the wall and Eo — 6 while in (b) 
the field is perpendicular to the wall and Eo = 12 is larger. 
The maximal density near the wall is similar in both panels 
but the interface is much wider in the perpendicular case due 
to favored head to tail configurations. 



field is obtained in the zero curvature limit for a capacitor 
consisting of concentric cylinders jl]. 

It is seen in Fig. [3] that for Eq = 4 the density pro- 
file shows a clear condensate in the parallel case. How- 
ever, for the same value of E in the perpendicular case, 
(p(z)} exhibits only a slight increase in the fluid density 
near the wall. Here, the field introduces a competition 
between alignment of dipoles parallel to the field, giving 
rise to the favored head-to-tail configurations, and the 
creation of an interface parallel to the field which dis- 
rupts these configurations [33) . This is in accord with 
the mean field description in which the typical field re- 
quired to induce condensation is an order of magnitude 
larger in the cylindrical capacitor compared to the wedge 
capacitor [4]. Only upon further increase of Eq to large 
values of E > 8 a significant increase in the density oc- 
curs. This is accompanied by large oscillations of the 
density close to the wall in which the distance between 
peaks is a. This is typical when fields perpendicular to 
the confining walls are applied to a high density dipolar 
fluid [33]. The oscillatory domain is followed by an in- 
terfacial region of width w 10cr, which is large compared 
to the thin interface of width k, 4er for parallel fields. 
Simulation snapshots of a small segment of the system, 
shown in Fig. [5j illustrate how the orientational order 
leads to a wider interface in the perpendicular case. We 
assume that due to the large interfacial energetic penalty 
in perpendicular fields one must simulate larger systems 
in order to observe clearly field induced condensation in 
this case. 

A hallmark of the first order transition in nonuni- 
form fields is a discontinuity in the surface density [5J. 
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FIG. 6. Maximal density as a function of the nonuniform field 
magnitude. Ao — 0.1 in all curves while the average density 
or number of particles is varied. 



The corresponding quantity in the simulation p max — 
max((p(z))) naturally occurs close to the wall where the 
field is large. In Fig. [6jwe plot p max as a function of the 
field magnitude. We find that p m ax{Eo) has a sigmoid 
like shape, similar to the mean field theory (see Fig. 10 
in Ref. [6j). However, since the simulated system is finite 
Pmax{Eo) changes continuously. Nonetheless, p m ax(E>o) 
grows more rapidly when the number of particles is in- 
creased from N = 8192 (squares) to N = 16384 (circles), 
suggesting that in the thermodynamic limit a first order 
transition is realized. 

Fig. |6j also shows that increasing the average density 
to p = 0.07 (diamonds) results in a larger condensate 
density. Although p = 0.07 is inside the binodal for the 
bulk system, this curve shows that one can utilize the 
nonuniform field to modify the density profile in a dilute 
enough finite system, such as a colloidal suspension. 

Results for the average dipole moment in the direction 
of the field (p x (z)} are shown in the solid curves of Fig. 
[71 These results are contrasted with the Debye theory 
[35] for an ideal gas shown in the dash-dot curves of Fig. 
[7] In the Debye theory: 



{»x{z)) = pL(a) 



(8) 



where L(a) = coth(a) — 1/a is the Langevin function and 
a = jit |-E(z)| /T. The simulation results agree with the 
Debye theory in the dilute vapor region but deviate to 
higher values in the dense liquid region. The deviation 
stems from the oversimplified treatment of the dipoles 
orientation correlation in the Debye theory |36j as well as 
the unaccounted effect of the short range LJ interaction 
on the orientational correlation. 

Further insight to the effect of the nonuniform field 
is gained by examining the polarization of the system. 
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FIG. 7. Solid curves: profiles of the scaled x component of the 
dipole moment fi x / for several nonuniform field strengths, 
dashed curves: Debye theory prediction (Eq. (|8j). 



In order to compare results for a uniform field with 
those of a nonuniform field we plot in the latter case the 
polarization as a function the average field 

1 f D 

Eq = -J^ E ( z ) dz ( 10 ) 

The solid curve in Fig. [8] shows that the polarization 
for the averaged nonuniform field is similar to that of 
the uniform field up to Eq sa 1. This value corresponds 
to Eq 2.5 which in Fig. [4] is where the fluid density 
near the wall starts to increase. For Eq > 1 the polar- 
ization rapidly increases as the fluid condensates until it 
saturates at large fields. Hence, the field induced con- 
densation can be utilized to amplify the electric response 
of a dilute dipolar system that will otherwise follow the 
Langevin type response. 
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FIG. 8. Scaled polarization P/P sa t as a function of the aver- 
age field E. Bulk and slab systems in a uniform field exhibit 
a Langevin type polarization while in a nonuniform field the 
polarization is larger starting at Bo K 1- 



follows from Eq. ^ that for E x = const.: 

(P) = P sat L{a) (9) 

where P sa t = Nfi is the saturation polarization of the 
system. We compare the polarization for a uniform field 
in the bulk and in the slab in Fig. [8] Simulation results 
in both cases are almost identical and agreement with the 
Debye theory is very good. This indicates that the bulk 
and slab system's response to a uniform field is essentially 
the same here because we consider a large enough slab, 
where the wall effects are negligible. 



IV. CONCLUSIONS 

We studied the effect of a nonuniform field on a Stock- 
mayer fluid via molecular dynamics simulations. We find 
that a homogeneous vapor phase in the canonical ensem- 
ble, unperturbed by a uniform field, undergoes a signifi- 
cant structural change in a nonuniform field of the same 
magnitude. This results in a sharp interface separating 
a liquid like region in the strong field region and a dilute 
vapor where the field is weaker. We attribute this change 
to the nonuniform field pulling the dipoles towards the 
strong field region combined with the attractive short 
range part of the potential. 

Our results indicate that a nonuniform field can be 
used to quite sensitively control the density profile and 
hence the fluid properties also in small closed systems. 
The mechanism we describe should be applicable for a 
broad class of one-component systems, including molecu- 
lar fluids and colloidal suspensions. In fact, a nonuniform 
field should promote phase separation in any dipolar sys- 
tem with an inherent bistable nature [12j [13] . We 
therefore believe that the study of fluids in nonuniform 
fields merits further experimental and theoretical atten- 
tion. 
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